Dataset downloaded from mgandal’s github repository.
# Load csvs
datExpr = read.csv('./../Data/RNAseq_ASD_datExpr.csv', row.names=1)
datMeta = read.csv('./../Data/RNAseq_ASD_datMeta.csv')
# 1. Group brain regions by lobes
# 2. Remove '/' from Batch variable: (It is recommended (but not required) to use only letters, numbers,
# and delimiters '_' or '.', in levels of factors as these are safe characters for column names in R
# 3. Transform Diagnosis into a factor variable
datMeta = datMeta %>% mutate(Brain_Region = as.factor(Region)) %>%
mutate(Brain_lobe = ifelse(Brain_Region %in% c('BA4_6', 'BA9', 'BA24', 'BA44_45'), 'Frontal',
ifelse(Brain_Region %in% c('BA3_1_2_5', 'BA7'), 'Parietal',
ifelse(Brain_Region %in% c('BA38', 'BA39_40', 'BA20_37', 'BA41_42_22'), 'Temporal',
'Occipital')))) %>%
mutate(Batch = as.factor(gsub('/', '.', RNAExtractionBatch)),
Diagnosis = factor(Diagnosis_, levels=c('CTL','ASD'))) %>%
dplyr::select(-Diagnosis_)
# Filter to only keep Frontal Lobe samples
datMeta = datMeta %>% filter(Brain_lobe=='Occipital')
datExpr = datExpr %>% dplyr::select(which(colnames(datExpr) %in% paste0('X',gsub('-','_',datMeta$X))))
# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>%
mutate('ID'=as.character(ensembl_gene_id)) %>%
dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
mutate('Neuronal'=1)
# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]
# NCBI biotype annotation
NCBI_biotype = read.csv('./../../../NCBI/Data/gene_biotype_info.csv') %>%
rename(Ensembl_gene_identifier='ensembl_gene_id', type_of_gene='gene_biotype', Symbol='hgnc_symbol') %>%
mutate(gene_biotype = ifelse(gene_biotype=='protein-coding','protein_coding',gene_biotype))
rm(GO_annotations)
Data description taken from the dataset’s synapse entry: RNAseq data was generated from 88 postmortem cortex brain samples from subjects with ASD (53 samples from 24 subjects) and non-psychiatric controls (35 samples from 17 subjects), across four cortical regions encompassing all major cortical lobes – frontal, temporal, parietal, and occipital. Brain samples were obtained from the Harvard Brain Bank as part of the Autism Tissue Project (ATP).
cat(paste0('Dataset includes ', nrow(datExpr), ' genes from ', ncol(datExpr), ' samples belonging to ',
length(unique(datMeta$Subject_ID)), ' different subjects.'))
## Dataset includes 63682 genes from 23 samples belonging to 23 different subjects.
Counts distribution: More than half of the counts are zero and most of the counts are relatively low, but there are some very high outliers
count_distr = summary(melt(datExpr))[,2]
for(i in 1:6){
print(count_distr[i])
}
##
## "Min. : 0 "
##
## "1st Qu.: 0 "
##
## "Median : 0 "
##
## "Mean : 565 "
##
## "3rd Qu.: 26 "
##
## "Max. :25859990 "
rm(count_distr, i)
Diagnosis distribution by Sample: There are more ASD samples than controls
table_info = datMeta %>% apply_labels(Diagnosis = 'Diagnosis', Brain_lobe = 'Brain Lobe', Batch = 'Batch', Sex = 'Gender')
cro(table_info$Diagnosis)
| #Total | |
|---|---|
| Diagnosis | |
| CTL | 8 |
| ASD | 15 |
| #Total cases | 23 |
Sex distribution: There are many more Male samples than Female ones
cro(table_info$Sex)
| #Total | |
|---|---|
| Gender | |
| F | 5 |
| M | 18 |
| #Total cases | 23 |
Diagnosis and Gender seem to be relatively balanced
cro(table_info$Diagnosis, list(table_info$Sex, total()))
| Gender | #Total | |||
|---|---|---|---|---|
| F | M | |||
| Diagnosis | ||||
| CTL | 1 | 7 | 8 | |
| ASD | 4 | 11 | 15 | |
| #Total cases | 5 | 18 | 23 | |
Age distribution: Subjects between 2 and 56 years old with a mean of 27
summary(datMeta$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.00 16.50 24.00 26.61 35.00 54.00
I was originally running this with the feb2014 version of BioMart because that’s the one that Gandal used (and it finds all of the Ensembl IDs, which other versions don’t), but it has some outdated biotype annotations, to fix this I’ll obtain all the information except the biotype label from BioMart in the same way as it had been done before, and then I’ll add the most current biotype label using information from NCBI’s website and then from BioMart in the following way:
1.Use BioMart to run a query with the original feb2014 version using the Ensembl IDs as keys to obtain all the information except the biotype labels
2.1 Use the NCBI annotations downloaded from NCBI’s website and processed in NCBI/RMarkdowns/20_02_07_clean_data.html (there is information for only 26K genes, so some genes will remain unlabelled)
2.2 Use the current version (jan2020) to obtain the biotype annotations using the Ensembl ID as keys (some genes don’t return a match)
2.3 For the genes that didn’t return a match, use the current version (jan2020) to obtain the biotype annotations using the gene name as keys (17 genes return multiple labels)
2.4 For the genes that returned multiple labels, use the feb2014 version with the Ensembl IDs as keys
labels_source = data.frame(data.frame('source' = c('NCBI', 'BioMart2020_byID', 'BioMart2020_byGene', 'BioMart2014'),
'n_matches' = rep(0,4)))
########################################################################################
# 1. Query archive version
getinfo = c('ensembl_gene_id','external_gene_id','chromosome_name','start_position',
'end_position','strand')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', host='feb2014.archive.ensembl.org')
datGenes = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=rownames(datExpr), mart=mart) %>%
rename(external_gene_id = 'hgnc_symbol')
## Cache found
datGenes$length = datGenes$end_position-datGenes$start_position
cat(paste0('1. ', sum(is.na(datGenes$start_position)), '/', nrow(datGenes),
' Ensembl IDs weren\'t found in the feb2014 version of BioMart'))
## 1. 0/63677 Ensembl IDs weren't found in the feb2014 version of BioMart
########################################################################################
########################################################################################
# 2. Get Biotype Labels
cat('2. Add biotype information')
## 2. Add biotype information
########################################################################################
# 2.1 Add NCBI annotations
datGenes = datGenes %>% left_join(NCBI_biotype, by=c('ensembl_gene_id','hgnc_symbol'))
cat(paste0('2.1 ' , sum(is.na(datGenes$gene_biotype)), '/', nrow(datGenes),
' Ensembl IDs weren\'t found in the NCBI database'))
## 2.1 42904/63677 Ensembl IDs weren't found in the NCBI database
labels_source$n_matches[1] = sum(!is.na(datGenes$gene_biotype))
########################################################################################
# 2.2 Query current BioMart version for gene_biotype using Ensembl ID as key
getinfo = c('ensembl_gene_id','gene_biotype')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', host='jan2020.archive.ensembl.org')
datGenes_biotype = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), mart=mart,
values=datGenes$ensembl_gene_id[is.na(datGenes$gene_biotype)])
## Cache found
cat(paste0('2.2 ' , sum(is.na(datGenes$gene_biotype))-nrow(datGenes_biotype), '/', sum(is.na(datGenes$gene_biotype)),
' Ensembl IDs weren\'t found in the jan2020 version of BioMart when querying by Ensembl ID'))
## 2.2 9099/42904 Ensembl IDs weren't found in the jan2020 version of BioMart when querying by Ensembl ID
# Add new gene_biotype info to datGenes
datGenes = datGenes %>% left_join(datGenes_biotype, by='ensembl_gene_id') %>%
mutate(gene_biotype = coalesce(as.character(gene_biotype.x), gene_biotype.y)) %>%
dplyr::select(-gene_biotype.x, -gene_biotype.y)
labels_source$n_matches[2] = sum(!is.na(datGenes$gene_biotype)) - labels_source$n_matches[1]
########################################################################################
# 3. Query current BioMart version for gene_biotype using gene symbol as key
missing_genes = unique(datGenes$hgnc_symbol[is.na(datGenes$gene_biotype)])
getinfo = c('hgnc_symbol','gene_biotype')
datGenes_biotype_by_gene = getBM(attributes=getinfo, filters=c('hgnc_symbol'), mart=mart,
values=missing_genes)
## Cache found
cat(paste0('2.3 ', length(missing_genes)-length(unique(datGenes_biotype_by_gene$hgnc_symbol)),'/',length(missing_genes),
' genes weren\'t found in the current BioMart version when querying by gene name'))
## 2.3 5712/7866 genes weren't found in the current BioMart version when querying by gene name
dups = unique(datGenes_biotype_by_gene$hgnc_symbol[duplicated(datGenes_biotype_by_gene$hgnc_symbol)])
cat(paste0(' ', length(dups), ' genes returned multiple labels (these won\'t be added)'))
## 17 genes returned multiple labels (these won't be added)
# Update information
datGenes_biotype_by_gene = datGenes_biotype_by_gene %>% filter(!hgnc_symbol %in% dups)
datGenes = datGenes %>% left_join(datGenes_biotype_by_gene, by='hgnc_symbol') %>%
mutate(gene_biotype = coalesce(gene_biotype.x, gene_biotype.y)) %>%
dplyr::select(-gene_biotype.x, -gene_biotype.y)
labels_source$n_matches[3] = sum(!is.na(datGenes$gene_biotype)) - sum(labels_source$n_matches)
########################################################################################
# 4. Query feb2014 BioMart version for the missing biotypes
missing_ensembl_ids = unique(datGenes$ensembl_gene_id[is.na(datGenes$gene_biotype)])
getinfo = c('ensembl_gene_id','gene_biotype')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', host='feb2014.archive.ensembl.org')
datGenes_biotype_archive = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=missing_ensembl_ids, mart=mart)
## Cache found
cat(paste0('2.4 ', length(missing_ensembl_ids)-nrow(datGenes_biotype_archive),'/',length(missing_ensembl_ids),
' genes weren\'t found in the feb2014 BioMart version when querying by Ensembl ID'))
## 2.4 0/6648 genes weren't found in the feb2014 BioMart version when querying by Ensembl ID
# Update information
datGenes = datGenes %>% left_join(datGenes_biotype_archive, by='ensembl_gene_id') %>%
mutate(gene_biotype = coalesce(gene_biotype.x, gene_biotype.y)) %>%
dplyr::select(-gene_biotype.x, -gene_biotype.y)
labels_source$n_matches[4] = sum(!is.na(datGenes$gene_biotype)) - sum(labels_source$n_matches)
########################################################################################
# Plot results
labels_source = labels_source %>% mutate(x = 1, percentage = round(100*n_matches/sum(n_matches),1))
p = labels_source %>% ggplot(aes(x, percentage, fill=source)) + geom_bar(position = 'stack', stat = 'identity') +
theme_minimal() + coord_flip() + theme(legend.position='bottom', axis.title.y=element_blank(),
axis.text.y=element_blank(), axis.ticks.y=element_blank())
ggplotly(p + theme(legend.position='none'))
as_ggplot(get_legend(p))
########################################################################################
# Reorder rows to match datExpr
datGenes = datGenes[match(rownames(datExpr), datGenes$ensembl_gene_id),]
rm(getinfo, mart, datGenes_biotype, datGenes_biotype_by_gene, datGenes_biotype_archive,
dups, missing_ensembl_ids, missing_genes, labels_source, p)
Checking how many SFARI genes are in the dataset
df = SFARI_genes %>% dplyr::select(-gene_biotype) %>% inner_join(datGenes, by=c('ID'='ensembl_gene_id'))
cat(paste0('Considering all genes, this dataset contains ', length(unique(df$`gene-symbol`)),
' of the ', length(unique(SFARI_genes$`gene-symbol`)), ' SFARI genes\n\n'))
## Considering all genes, this dataset contains 979 of the 980 SFARI genes
cat(paste0('The missing gene is ',
SFARI_genes$`gene-symbol`[! SFARI_genes$`gene-symbol` %in% df$`gene-symbol`],
' with a SFARI score of ',
SFARI_genes$`gene-score`[! SFARI_genes$`gene-symbol` %in% df$`gene-symbol`]))
## The missing gene is MIR137 with a SFARI score of 3
rm(df)
1. Filter entries that don’t correspond to genes
to_keep = !is.na(datGenes$length)
cat(paste0('Names of the rows removed: ', paste(rownames(datExpr)[!to_keep], collapse=', ')))
## Names of the rows removed: __no_feature, __ambiguous, __too_low_aQual, __not_aligned, __alignment_not_unique
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
rownames(datGenes) = datGenes$ensembl_gene_id
cat(paste0('Removed ', sum(!to_keep), ' \'genes\', ', sum(to_keep), ' remaining'))
## Removed 5 'genes', 63677 remaining
2. Filter genes that do not encode any protein
cat(paste0(sum(datGenes$gene_biotype=='protein_coding'), '/', nrow(datGenes), ' are protein coding genes' ))
## 22543/63677 are protein coding genes
sort(table(datGenes$gene_biotype), decreasing=TRUE)
##
## protein_coding lncRNA
## 22543 12167
## processed_pseudogene unprocessed_pseudogene
## 10117 2547
## 1 miRNA
## 2314 2276
## misc_RNA snRNA
## 2178 2043
## pseudogene snoRNA
## 1410 1202
## lincRNA transcribed_unprocessed_pseudogene
## 840 682
## rRNA_pseudogene transcribed_processed_pseudogene
## 500 441
## antisense 3
## 380 331
## 6 IG_V_pseudogene
## 314 254
## IG_V_gene TR_V_gene
## 179 146
## transcribed_unitary_pseudogene TR_J_gene
## 86 81
## unitary_pseudogene processed_transcript
## 74 72
## sense_intronic IG_D_gene
## 72 64
## rRNA TR_V_pseudogene
## 49 46
## sense_overlapping scaRNA
## 38 31
## polymorphic_pseudogene 7
## 28 25
## IG_J_gene IG_C_gene
## 24 23
## Mt_tRNA 4
## 22 17
## IG_C_pseudogene TEC
## 11 11
## TR_C_gene 3prime_overlapping_ncrna
## 8 6
## IG_J_pseudogene ribozyme
## 6 5
## TR_J_pseudogene TR_D_gene
## 5 3
## Mt_rRNA 8
## 2 1
## translated_processed_pseudogene translated_unprocessed_pseudogene
## 1 1
## vaultRNA
## 1
Most of the genes with low expression levels are not protein-coding
plot_data = data.frame('ID' = rownames(datExpr), 'MeanExpr' = apply(datExpr, 1, mean),
'ProteinCoding' = datGenes$gene_biotype=='protein_coding')
ggplotly(plot_data %>% ggplot(aes(log2(MeanExpr+1), fill=ProteinCoding, color=ProteinCoding)) + geom_density(alpha=0.5) +
theme_minimal())
rm(plot_data)
We lose 3 genes with a SFARI score
Note: The gene name for Ensembl ID ENSG00000187951 is wrong, it should be AC091057.1 instead of ARHGAP11B, but the biotype is right, so it would still be filtered out
df = SFARI_genes %>% dplyr::select(-gene_biotype) %>% inner_join(datGenes, by=c('ID'='ensembl_gene_id'))
cat(paste0('Filtering protein coding genes, we are left with ', length(unique(df$`gene-symbol`[df$gene_biotype=='protein_coding'])),
' SFARI genes'))
## Filtering protein coding genes, we are left with 976 SFARI genes
kable(df %>% filter(! `gene-symbol` %in% df$`gene-symbol`[df$gene_biotype=='protein_coding']) %>%
dplyr::select(ID, `gene-symbol`, `gene-score`, gene_biotype, syndromic, `number-of-reports`), caption='Lost Genes')
| ID | gene-symbol | gene-score | gene_biotype | syndromic | number-of-reports |
|---|---|---|---|---|---|
| ENSG00000187951 | ARHGAP11B | 4 | lncRNA | 0 | 2 |
| ENSG00000251593 | MSNP1AS | 2 | processed_pseudogene | 0 | 12 |
| ENSG00000197558 | SSPO | 4 | transcribed_unitary_pseudogene | 0 | 3 |
rm(df)
if(!all(rownames(datExpr)==rownames(datGenes))) cat('!!! gene rownames do not match!!!')
to_keep = datGenes$gene_biotype=='protein_coding'
datExpr = datExpr %>% filter(to_keep)
datGenes = datGenes %>% filter(to_keep)
rownames(datExpr) = datGenes$ensembl_gene_id
rownames(datGenes) = datGenes$ensembl_gene_id
cat(paste0(length(unique(SFARI_genes$`gene-symbol`[SFARI_genes$ID %in% rownames(datExpr)])), ' SFARI genes remaining'))
## 976 SFARI genes remaining
cat(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## Removed 41134 genes, 22543 remaining
3. Filter genes with low expression levels
\(\qquad\) 3.1 Remove genes with zero expression in all of the samples
to_keep = rowSums(datExpr)>0
cat(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## Removed 3830 genes, 18713 remaining
df = data.frame('rowSums' = rowSums(datExpr), 'ensembl_gene_id' = rownames(datExpr)) %>%
right_join(SFARI_genes, by='ensembl_gene_id') %>% filter(rowSums==0 & `gene-score` %in% c(1,2,3)) %>%
arrange(`gene-score`) %>% dplyr::select(-ensembl_gene_id) %>% filter(!duplicated(`gene-symbol`))
kable(df %>% dplyr::select(ID, `gene-symbol`, `gene-score`, gene_biotype, syndromic, `number-of-reports`),
caption='Lost Genes with Top Scores')
| ID | gene-symbol | gene-score | gene_biotype | syndromic | number-of-reports |
|---|---|---|---|---|---|
| ENSG00000267910 | KMT2A | 1 | protein_coding | 1 | 20 |
| ENSG00000227460 | SYNGAP1 | 1 | protein_coding | 1 | 50 |
| ENSG00000265594 | CEP41 | 2 | protein_coding | 0 | 5 |
| ENSG00000272883 | CNOT3 | 2 | protein_coding | 1 | 5 |
| ENSG00000259938 | CUX1 | 2 | protein_coding | 0 | 6 |
| ENSG00000271019 | MBOAT7 | 2 | protein_coding | 1 | 3 |
| ENSG00000268563 | MECP2 | 2 | protein_coding | 1 | 75 |
| ENSG00000262024 | TCF20 | 2 | protein_coding | 1 | 10 |
| ENSG00000266334 | APH1A | 3 | protein_coding | 0 | 2 |
| ENSG00000260508 | GGNBP2 | 3 | protein_coding | 0 | 2 |
| ENSG00000272706 | GRIK5 | 3 | protein_coding | 0 | 8 |
| ENSG00000269816 | KDM5C | 3 | protein_coding | 0 | 22 |
| ENSG00000235718 | MFRP | 3 | protein_coding | 0 | 6 |
| ENSG00000267946 | TMLHE | 3 | protein_coding | 0 | 5 |
| ENSG00000260150 | ZMYND11 | 3 | protein_coding | 0 | 10 |
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
cat(paste0(length(unique(SFARI_genes$`gene-symbol`[SFARI_genes$ID %in% rownames(datExpr)])), ' SFARI genes remaining'))
## 964 SFARI genes remaining
rm(df)
\(\qquad\) 2.2 Removing genes with a high percentage of zeros
Choosing the threshold:
Criteria for selecting the percentage of zeros threshold: The minimum value in which the preprocessed data is relatively homoscedastic (we’re trying to get rid of the group of genes with very low mean and SD that make the cloud of points look like a comic book speech bubble)
On the plot I’m using the “dual” of the maximum percentage of zeros, the minimum percentage of non zeros so the visualisation is more intuitive
57% seems to be a good threshold for the minimum percentage of non zeros, so 43% will be the maximum percentage of zeros allowed in a row
The Mean vs SD plot doesn’t show all of the genes, a random sample was selected for the genes with higher level of expression so the visualisation wouldn’t be as heavy (and since we care about the genes with the lowest levels of expression, we aren’t losing important information)
datMeta_original = datMeta
datExpr_original = datExpr
datGenes_original = datGenes
# Return to original variables
datExpr = datExpr_original
datGenes = datGenes_original
datMeta = datMeta_original
rm(datExpr_original, datGenes_original, datMeta_original, datExpr_vst, datGenes_vst, datMeta_vst)
Filtering
# Minimum percentage of non-zero entries allowed per gene
threshold = 57
plot_data = data.frame('id'=rownames(datExpr), 'non_zero_percentage' = apply(datExpr, 1, function(x) 100*mean(x>0)))
ggplotly(plot_data %>% ggplot(aes(x=non_zero_percentage)) + geom_density(color='#0099cc', fill='#0099cc', alpha=0.3) +
geom_vline(xintercept=threshold, color='gray') + #scale_x_log10() +
ggtitle('Percentage of non-zero entries distribution') + theme_minimal())
to_keep = apply(datExpr, 1, function(x) 100*mean(x>0)) >= threshold
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
cat(paste0(length(unique(SFARI_genes$`gene-symbol`[SFARI_genes$ID %in% rownames(datExpr)])), ' SFARI genes remaining'))
## 931 SFARI genes remaining
cat(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## Removed 2223 genes, 16490 remaining
rm(threshold, plot_data, to_keep)
3. Filter outlier samples
\(\qquad\) 3.1 Gandal filters samples belonging to subject AN03345 without giving an explanation. Since it could have some technical problems, I remove them as well
to_keep = (datMeta$Subject_ID != 'AN03345')
datMeta = datMeta[to_keep,]
datExpr = datExpr[,to_keep]
cat(paste0('Removed ', sum(!to_keep), ' samples, ', sum(to_keep), ' remaining'))
## Removed 0 samples, 23 remaining
\(\qquad\) 3.2 Filter out outliers: Using node connectivity as a distance measure, normalising it and filtering out genes farther away than 2 standard deviations from the left (lower connectivity than average, not higher)
Gandal uses the formula \(s_{ij}=\frac{1+bw(i,j)}{2}\) to convert all the weights to positive values, but I used \(s_{ij}=|bw(i,j)|\) instead because I think it makes more sense. In the end it doesn’t matter because they select as outliers the same six samples
Two outlier samples were found, both from the ASD group
absadj = datExpr %>% bicor %>% abs
## alpha: 1.000000
netsummary = fundamentalNetworkConcepts(absadj)
ku = netsummary$Connectivity
z.ku = (ku-mean(ku))/sqrt(var(ku))
plot_data = data.frame('sample'=1:length(z.ku), 'distance'=z.ku, 'Sample_ID'=datMeta$Sample_ID,
'Subject_ID'=datMeta$Subject_ID, 'Extraction_Batch'=datMeta$RNAExtractionBatch,
'Brain_Lobe'=datMeta$Brain_lobe, 'Sex'=datMeta$Sex, 'Age'=datMeta$Age,
'Diagnosis'=datMeta$Diagnosis, 'PMI'=datMeta$PMI)
selectable_scatter_plot(plot_data, plot_data[,-c(1,2)])
cat(paste0('Outlier sample(s): ', paste(as.character(plot_data$Sample_ID[plot_data$distance< -2]), collapse=', ')))
## Outlier sample(s): AN17254_BA17, AN02987_BA17
to_keep = z.ku > -2
datMeta = datMeta[to_keep,]
datExpr = datExpr[,to_keep]
cat(paste0('Removed ', sum(!to_keep), ' samples, ', sum(to_keep), ' remaining'))
## Removed 2 samples, 21 remaining
rm(absadj, netsummary, ku, z.ku, plot_data, to_keep)
cat(paste0('After filtering, the dataset consists of ', nrow(datExpr), ' genes and ', ncol(datExpr), ' samples'))
## After filtering, the dataset consists of 16490 genes and 21 samples
According to Tackling the widespread and critical impact of batch effects in high-throughput data, technical artifacts can be an important source of variability in the data, so batch correction should be part of the standard preprocessing pipeline of gene expression data.
They say Processing group and Date of the experiment are good batch surrogates, so I’m going to see if they affect the data in any clear way to use them as surrogates.
All the information we have is the Brain Bank, and although all the samples were obtained from the Autism Tissue Project, we don’t have any more specific information about who preprocessed each sample
table(datMeta$Brain_Bank)
##
## ATP
## 21
There are two different dates when the data was procesed
table(datMeta$RNAExtractionBatch)
##
## 10/10/2014 6/20/2014
## 14 7
On the whole dataset there wasn’t a correlation betwen batch surrogante and diagnosis, but now that we are filtering for Occipital Lobe samples, all but one of the control samples were processed on the first batch
table(datMeta$RNAExtractionBatch, datMeta$Diagnosis)
##
## CTL ASD
## 10/10/2014 7 7
## 6/20/2014 1 6
Samples don’t seem to cluster together that strongly for each batch and I can’t see a confounding effect between batch and diagnosis, but perhaps it’s because there aren’t enough observations
h_clusts = datExpr %>% t %>% dist %>% hclust %>% as.dendrogram
create_viridis_dict = function(){
min_age = datMeta$Age %>% min
max_age = datMeta$Age %>% max
viridis_age_cols = viridis(max_age - min_age + 1)
names(viridis_age_cols) = seq(min_age, max_age)
return(viridis_age_cols)
}
viridis_age_cols = create_viridis_dict()
dend_meta = datMeta[match(substring(labels(h_clusts),2), datMeta$Dissected_Sample_ID),] %>%
mutate('Batch' = ifelse(RNAExtractionBatch=='10/10/2014', '#F8766D', '#00BFC4'),
'Diagnosis' = ifelse(Diagnosis=='CTL','#008080','#86b300'), # Blue control, Green ASD
'Sex' = ifelse(Sex=='F','#ff6666','#008ae6'), # Pink Female, Blue Male
'Age' = viridis_age_cols[as.character(Age)]) %>% # Purple: young, Yellow: old
dplyr::select(Age, Sex, Diagnosis, Batch)
h_clusts %>% as.dendrogram %>% dendextend::set('labels', rep('', nrow(datMeta))) %>%
dendextend::set('branches_k_color', k=5) %>% plot
colored_bars(colors=dend_meta)
rm(h_clusts, dend_meta, create_viridis_dict, viridis_age_cols)
Comparing the mean expression of each sample by batch, this time, the difference between batches doesn’t seem that big
plot_data_b1 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='10/10/2014']), 'Batch'='10/10/2014')
plot_data_b2 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='6/20/2014']), 'Batch'='6/20/2014')
plot_data = rbind(plot_data_b1, plot_data_b2)
mu = plot_data %>% group_by(Batch) %>% dplyr::summarise(BatchMean=mean(Mean))
ggplotly(plot_data %>% ggplot(aes(x=Mean, color=Batch, fill=Batch)) + geom_density(alpha=0.3) +
geom_vline(data=mu, aes(xintercept=BatchMean, color=Batch), linetype='dashed') +
ggtitle('Mean expression by sample grouped by Batch') + scale_x_log10() + theme_minimal())
rm(plot_data_b1, plot_data_b2, plot_data, mu)
Following the pipeline from Surrogate variable analysis: hidden batch effects where sva is used with DESeq2.
Create a DeseqDataSet object, estimate the library size correction and save the normalized counts matrix
counts = datExpr %>% as.matrix
rowRanges = GRanges(datGenes$chromosome_name,
IRanges(datGenes$start_position, width=datGenes$length),
strand=datGenes$strand,
feature_id=datGenes$ensembl_gene_id)
se = SummarizedExperiment(assays=SimpleList(counts=counts), rowRanges=rowRanges, colData=datMeta)
dds = DESeqDataSet(se, design = ~ Diagnosis)
dds = estimateSizeFactors(dds)
norm.cts = counts(dds, normalized=TRUE)
Provide the normalized counts and two model matrices to SVA. The first matrix uses the biological condition, and the second model matrix is the null model.
mod = model.matrix(~ Diagnosis, colData(dds))
mod0 = model.matrix(~ 1, colData(dds))
sva_fit = svaseq(norm.cts, mod=mod, mod0=mod0)
## Number of significant surrogate variables is: 4
## Iteration (out of 5 ):1 2 3 4 5
rm(mod, mod0, norm.cts)
Found 4 surrogate variables, since there is no direct way to select which ones to pick Bioconductor answer, kept all of them.
Include SV estimations to datMeta information
sv_data = sva_fit$sv %>% data.frame
colnames(sv_data) = paste0('SV',1:ncol(sv_data))
datMeta_sva = cbind(datMeta, sv_data)
rm(sv_data, sva_fit)
In conclusion: Date of extraction works as a surrogate for batch effect and the sva package found other 5 variables that could work as surrogates which are now included in datMeta and should be included in the DEA.
Using DESeq2 package to perform normalisation. I chose this package over limma because limma uses the log transformed data as input instead of the raw counts and I have discovered that in this dataset, this transformation affects genes differently depending on their mean expression level, and genes with a high SFARI score are specially affected by this.
plot_data = data.frame('ID'=rownames(datExpr), 'Mean'=rowMeans(datExpr), 'SD'=apply(datExpr,1,sd))
plot_data %>% ggplot(aes(Mean, SD)) + geom_point(color='#0099cc', alpha=0.1) + geom_abline(color='gray') +
scale_x_log10() + scale_y_log10() + theme_minimal()
rm(plot_data)
Using vst instead of rlog to perform normalisation. Bioconductor question explaining differences between methods. Chose vst because a) it is much faster than rlog (it is recommended to use vst for samples larger than 50), and b) Michael Love (author of DESEq2) recommends using it over rlog
Including a log fold change threshold of 0 in the results formula \(H_0:lfc=0\) because setting any other log fold change seems arbitrary and we risk losing genes with a significant differential expression for genes with a higher fold change, but not necessarily as significant.
counts = datExpr %>% as.matrix
rowRanges = GRanges(datGenes$chromosome_name,
IRanges(datGenes$start_position, width=datGenes$length),
strand=datGenes$strand,
feature_id=datGenes$ensembl_gene_id)
se = SummarizedExperiment(assays=SimpleList(counts=counts), rowRanges=rowRanges, colData=datMeta_sva)
dds = DESeqDataSet(se, design = ~ Batch + SV1 + SV2 + SV3 + SV4 + Diagnosis)
# Perform DEA
dds = DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DE_info = results(dds, lfcThreshold=0, altHypothesis='greaterAbs')
# Perform vst
vsd = vst(dds)
datExpr_vst = assay(vsd)
datMeta_vst = colData(vsd)
datGenes_vst = rowRanges(vsd)
rm(counts, rowRanges, se, vsd)
Using the plotting function DESEq2’s manual proposes to study vst’s output it looks like the data could be homoscedastic
meanSdPlot(datExpr_vst, plot=FALSE)$gg + theme_minimal()
When plotting point by point it seems like the group of genes with the lowest values behave differently to the rest
plot_data = data.frame('ID'=rownames(datExpr_vst), 'Mean'=rowMeans(datExpr_vst), 'SD'=apply(datExpr_vst,1,sd))
plot_data %>% ggplot(aes(Mean, SD)) + geom_point(color='#0099cc', alpha=0.2) +
scale_x_log10() + scale_y_log10() + theme_minimal()
rm(plot_data)
*Could have done this since before
save(datExpr, datMeta, datGenes, file='./../Data/filtered_raw_data.RData')
#load('./../Data/filtered_raw_data.RData')
Rename normalised datasets to continue working with these
datExpr = datExpr_vst
datMeta = datMeta_vst %>% data.frame
datGenes = datGenes_vst
cat(paste0(length(unique(SFARI_genes$`gene-symbol`[SFARI_genes$ID %in% rownames(datExpr)])), ' SFARI genes remaining'))
## 931 SFARI genes remaining
cat(paste0('After filtering, the dataset consists of ', nrow(datExpr), ' genes and ', ncol(datExpr), ' samples'))
## After filtering, the dataset consists of 16490 genes and 21 samples
rm(datExpr_vst, datMeta_vst, datGenes_vst, datMeta_sva)
By including the surrogate variables in the DESeq formula we only modelled the batch effects into the DEA, but we didn’t actually correct them from the data, for that we need to use ComBat (or other equivalent package) in the already normalised data
In some places they say you shouldn’t correct these effects on the data because you risk losing biological variation, in others they say you should because they introduce noise to the data. The only thing everyone agrees on is that you shouldn’t remove them before performing DEA but instead include them in the model.
Based on the conclusions from Practical impacts of genomic data “cleaning” on biological discovery using surrogate variable analysis it seems like it may be a good idea to remove the batch effects from the data and not only from the DE analysis:
Using SVA, ComBat or related tools can increase the power to identify specific signals in complex genomic datasets (they found “greatly sharpened global and gene-specific differential expression across treatment groups”)
But caution should be exercised to avoid removing biological signal of interest
We must be precise and deliberate in the design and analysis of experiments and the resulting data, and also mindful of the limitations we impose with our own perspective
Open data exploration is not possible after such supervised “cleaning”, because effects beyond those stipulated by the researcher may have been removed
# Taken from https://www.biostars.org/p/121489/#121500
correctDatExpr = function(datExpr, mod, svs) {
X = cbind(mod, svs)
Hat = solve(t(X) %*% X) %*% t(X)
beta = (Hat %*% t(datExpr))
rm(Hat)
gc()
P = ncol(mod)
return(datExpr - t(as.matrix(X[,-c(1:P)]) %*% beta[-c(1:P),]))
}
pca_samples_before = datExpr %>% t %>% prcomp
pca_genes_before = datExpr %>% prcomp
# Correct
mod = model.matrix(~ Diagnosis, colData(dds))
svs = datMeta %>% dplyr::select(SV1:SV4) %>% as.matrix
datExpr_corrected = correctDatExpr(as.matrix(datExpr), mod, svs)
pca_samples_after = datExpr_corrected %>% t %>% prcomp
pca_genes_after = datExpr_corrected %>% prcomp
rm(correctDatExpr)
Removing batch effects has a big impact in the distribution of the samples, unlike the Frontal Lobe samples, they manage to separate perfectly by the first principal component after the SVA correction
pca_samples_df = rbind(data.frame('ID'=colnames(datExpr), 'PC1'=pca_samples_before$x[,1],
'PC2'=pca_samples_before$x[,2], 'corrected'=0),
data.frame('ID'=colnames(datExpr), 'PC1'=pca_samples_after$x[,1],
'PC2'=pca_samples_after$x[,2], 'corrected'=1)) %>%
left_join(datMeta %>% mutate('ID'=rownames(datMeta)), by='ID')
ggplotly(pca_samples_df %>% ggplot(aes(PC1, PC2, color=Diagnosis)) + geom_point(aes(frame=corrected, id=ID), alpha=0.75) +
xlab(paste0('PC1 (corr=', round(cor(pca_samples_before$x[,1],pca_samples_after$x[,1]),2),
'). % Var explained: ', round(100*summary(pca_samples_before)$importance[2,1],1),' to ',
round(100*summary(pca_samples_after)$importance[2,1],1))) +
ylab(paste0('PC2 (corr=', round(cor(pca_samples_before$x[,2],pca_samples_after$x[,2]),2),
'). % Var explained: ', round(100*summary(pca_samples_before)$importance[2,2],1),' to ',
round(100*summary(pca_samples_after)$importance[2,2],1))) +
ggtitle('Samples') + theme_minimal())
rm(pca_samples_df)
It seems like the sva correction preserves the mean expression of the genes and erases almost everything else (although what little else remains is enough to characterise the two Diagnosis groups pretty well using only the first PC)
*Plot is done with only 10% of the genes so it’s not that heavy
pca_genes_df = rbind(data.frame('ID'=rownames(datExpr), 'PC1'=pca_genes_before$x[,1],
'PC2'=pca_genes_before$x[,2], 'corrected'=0, 'MeanExpr'=rowMeans(datExpr)),
data.frame('ID'=rownames(datExpr), 'PC1'=pca_genes_after$x[,1],
'PC2'=pca_genes_after$x[,2], 'corrected'=1, 'MeanExpr'=rowMeans(datExpr)))
keep_genes = rownames(datExpr) %>% sample(0.1*nrow(datExpr))
pca_genes_df = pca_genes_df %>% filter(ID %in% keep_genes)
ggplotly(pca_genes_df %>% ggplot(aes(PC1, PC2,color=MeanExpr)) + geom_point(alpha=0.3, aes(frame=corrected, id=ID)) +
xlab(paste0('PC1 (corr=', round(cor(pca_genes_before$x[,1],pca_genes_after$x[,1]),2),
'). % Var explained: ', round(100*summary(pca_genes_before)$importance[2,1],1),' to ',
round(100*summary(pca_genes_after)$importance[2,1],1))) +
ylab(paste0('PC2 (corr=', round(cor(pca_genes_before$x[,2],pca_genes_after$x[,2]),2),
'). % Var explained: ', round(100*summary(pca_genes_before)$importance[2,2],1),' to ',
round(100*summary(pca_genes_after)$importance[2,2],1))) +
scale_color_viridis() + ggtitle('Genes') + theme_minimal())
rm(pca_samples_before, pca_genes_before, mod, svs, pca_samples_after, pca_genes_after, pca_genes_df, keep_genes)
Everything looks good, so we’re keeping the corrected expression dataset
datExpr = datExpr_corrected
rm(datExpr_corrected)
Even after correcting the dataset for the surrogate variables found with sva, there is still a difference in mean expression by processing date (although this difference is relatively small)
Since the processing date and Diagnosis are correlated and because we know from the exploratory analysis of the whole dataset that ASD samples tend to have higher levels of expression than the controls, we cannot tell if the difference we are seeing is coming from the Diagnosis or from the Batch effect (because the curve with the lower level of expression corresponds to the batch with the majority of the CTL samples). Because of this, I’m not going to correct for processing date.
plot_data_b1 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='10/10/2014']), 'Batch'='10/10/2014')
plot_data_b2 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='6/20/2014']), 'Batch'='6/20/2014')
plot_data = rbind(plot_data_b1, plot_data_b2)
mu = plot_data %>% group_by(Batch) %>% dplyr::summarise(BatchMean=mean(Mean))
ggplotly(plot_data %>% ggplot(aes(x=Mean, color=Batch, fill=Batch)) + geom_density(alpha=0.3) +
geom_vline(data=mu, aes(xintercept=BatchMean, color=Batch), linetype='dashed') +
ggtitle('Mean expression by sample grouped by processing date') + scale_x_log10() + theme_minimal())
rm(plot_data_b1, plot_data_b2, plot_data, mu)
save(datExpr, datMeta, datGenes, DE_info, dds, file='./../Data/preprocessed_data.RData')
#load('./../Data/preprocessed_data.RData')
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] knitr_1.24 expss_0.10.1
## [3] dendextend_1.13.3 vsn_3.54.0
## [5] WGCNA_1.68 fastcluster_1.1.25
## [7] dynamicTreeCut_1.63-1 sva_3.34.0
## [9] genefilter_1.68.0 mgcv_1.8-28
## [11] nlme_3.1-139 DESeq2_1.26.0
## [13] SummarizedExperiment_1.16.1 DelayedArray_0.12.2
## [15] BiocParallel_1.20.1 matrixStats_0.55.0
## [17] Biobase_2.46.0 GenomicRanges_1.38.0
## [19] GenomeInfoDb_1.22.0 IRanges_2.20.2
## [21] S4Vectors_0.24.3 BiocGenerics_0.32.0
## [23] biomaRt_2.42.0 ggpubr_0.2.5
## [25] magrittr_1.5 ggExtra_0.9
## [27] GGally_1.4.0 gridExtra_2.3
## [29] viridis_0.5.1 viridisLite_0.3.0
## [31] RColorBrewer_1.1-2 plotlyutils_0.0.0.9000
## [33] plotly_4.9.2 glue_1.3.1
## [35] reshape2_1.4.3 forcats_0.4.0
## [37] stringr_1.4.0 dplyr_0.8.3
## [39] purrr_0.3.3 readr_1.3.1
## [41] tidyr_1.0.2 tibble_2.1.3
## [43] ggplot2_3.2.1 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 robust_0.4-18.2 RSQLite_2.2.0
## [4] AnnotationDbi_1.48.0 htmlwidgets_1.5.1 grid_3.6.0
## [7] munsell_0.5.0 codetools_0.2-16 preprocessCore_1.48.0
## [10] miniUI_0.1.1.1 withr_2.1.2 colorspace_1.4-1
## [13] highr_0.8 rstudioapi_0.10 robustbase_0.93-5
## [16] ggsignif_0.6.0 labeling_0.3 GenomeInfoDbData_1.2.2
## [19] bit64_0.9-7 farver_2.0.3 vctrs_0.2.2
## [22] generics_0.0.2 xfun_0.8 BiocFileCache_1.10.2
## [25] R6_2.4.1 doParallel_1.0.15 locfit_1.5-9.1
## [28] bitops_1.0-6 reshape_0.8.8 assertthat_0.2.1
## [31] promises_1.1.0 scales_1.1.0 nnet_7.3-12
## [34] gtable_0.3.0 Cairo_1.5-10 affy_1.64.0
## [37] rlang_0.4.4 splines_3.6.0 lazyeval_0.2.2
## [40] acepack_1.4.1 impute_1.60.0 hexbin_1.28.1
## [43] broom_0.5.4 checkmate_1.9.4 BiocManager_1.30.10
## [46] yaml_2.2.0 modelr_0.1.5 crosstalk_1.0.0
## [49] backports_1.1.5 httpuv_1.5.2 Hmisc_4.2-0
## [52] tools_3.6.0 affyio_1.56.0 ellipsis_0.3.0
## [55] Rcpp_1.0.3 plyr_1.8.5 base64enc_0.1-3
## [58] progress_1.2.2 zlibbioc_1.32.0 RCurl_1.95-4.12
## [61] prettyunits_1.0.2 rpart_4.1-15 openssl_1.4.1
## [64] cowplot_1.0.0 haven_2.2.0 cluster_2.0.8
## [67] fs_1.3.1 data.table_1.12.8 reprex_0.3.0
## [70] mvtnorm_1.0-11 hms_0.5.3 mime_0.9
## [73] evaluate_0.14 xtable_1.8-4 XML_3.99-0.3
## [76] readxl_1.3.1 compiler_3.6.0 crayon_1.3.4
## [79] htmltools_0.4.0 pcaPP_1.9-73 later_1.0.0
## [82] Formula_1.2-3 geneplotter_1.64.0 rrcov_1.4-7
## [85] lubridate_1.7.4 DBI_1.1.0 dbplyr_1.4.2
## [88] MASS_7.3-51.4 rappdirs_0.3.1 Matrix_1.2-17
## [91] cli_2.0.1 pkgconfig_2.0.3 fit.models_0.5-14
## [94] foreign_0.8-71 xml2_1.2.2 foreach_1.4.7
## [97] annotate_1.64.0 XVector_0.26.0 rvest_0.3.5
## [100] digest_0.6.24 rmarkdown_1.14 cellranger_1.1.0
## [103] htmlTable_1.13.1 curl_4.3 shiny_1.4.0
## [106] lifecycle_0.1.0 jsonlite_1.6 askpass_1.1
## [109] limma_3.42.2 fansi_0.4.1 pillar_1.4.3
## [112] lattice_0.20-38 fastmap_1.0.1 httr_1.4.1
## [115] DEoptimR_1.0-8 survival_2.44-1.1 GO.db_3.10.0
## [118] iterators_1.0.12 bit_1.1-15.2 stringi_1.4.6
## [121] blob_1.2.1 latticeExtra_0.6-28 memoise_1.1.0